1 Importación de datos

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
calidad_aire <- readRDS("data/02_estructura.RDS")
calidad_aire_orig <- readRDS("data/raw/calidad_aire.RDS")
estaciones <- readRDS("data/raw/coordenadas_estaciones.RDS")

El conjunto de datos calidad del aire contiene 128357 observaciones y 19 variables. En las secciones siguientes analizaremos los datos haciendo las operaciones de transformación y limpieza necesarias.

2 Estaciones de medición

Existen 24 estaciones distintas cuya situación se puede ver en el siguiente mapa

library(leaflet)
leaflet(data = estaciones) %>%
  addTiles() %>%
  addMarkers(lng = ~longitud, lat = ~latitud,
             label = ~as.character(estacion)
             )

Podemos estudiar la calidad de los datos en las distintas estaciones. Por ejemplo, puede que no todas las estaciones sean capaces de medir todas las magnitudes. En el siguiente gráfico representamos el porcentaje de mediciones ausentes (o no válidas) para cada magnitud y estación

porc_na <- calidad_aire %>% 
  group_by(estacion) %>% 
  summarise_at(vars(so2:nmhc), function(x) mean(is.na(x)))

porc_na %>% 
  pivot_longer(
    cols = so2:nmhc,
    names_to = "magnitud"
  ) %>% 
  ggplot(
    aes(x = magnitud,
        y = as.factor(estacion),
        fill = value)
  ) +
  geom_tile() +
  theme_minimal() +
  labs(
    title = "Porcentaje de mediciones ausentes para cada magnitud y estación",
    fill = "",
    x = "",
    y = "Estación",
    caption = "MSMK: taller de calidad del aire"
  ) +
  scale_fill_continuous(trans = "reverse", labels = scales::percent, breaks = c(1, 0.5, 0), limits = c(1, 0)) +
  theme(legend.position = "top")

En el gráfico anterior se puede ver que la calidad de la información de las magnitudes no, no2 y nox es muy alta. Sin embargo, para el resto de magnitudes no hay homogeneidad entre las estaciones. Además, en este proyecto buscamos predecir la calidad del aire para la ciudad de Madrid. Por lo tanto, la información deberá estar agregada de forma que tengamos un único valor por magnitud y día.

La pregunta es ¿cómo agregamos la información? La aproximación más habitual sería obtener el valor medio de cada magnitud en cada día. Ésta es la agregación que seguiremos en este proyecto. No obstante, la agregación es una cuestión importante en el diseño del objetivo del proyecto. Por ejemplo, el valor máximo podría ser otra forma de agregación, de forma que se pudiesen activar protocolos de contaminación en el caso de que en algún punto de la ciudad se sobrepasen ciertos límites y que, utilizando la media, podríamos diluir los valores extremos.

Agrupamos los datos con la media en el siguiente código:

calidad_aire <- calidad_aire %>% 
  group_by(fecha, anyo, mes, dia) %>% 
  summarise_at(vars(so2:nmhc), mean, na.rm = TRUE) %>% 
  ungroup()

Nos podemos hacer la pregunta de si la calidad de los datos ha ido evolucionando con el tiempo. A continuación representamos aquellos días para los que las magnitudes no están informadas aún habiendo agregado todas las estaciones

calidad_aire %>% 
  pivot_longer(
    cols = so2:nmhc,
    names_to = "magnitud"
  ) %>% 
  mutate(value = is.na(value)) %>% 
  filter(value == TRUE) %>% 
  ggplot(
    aes(x = magnitud,
        y = fecha,
        fill = as.character(as.numeric(value))
    )
  ) +
  geom_tile() +
  scale_fill_manual(values = c("1" = "black", "0" = "white")) +
  theme_minimal() +
  labs(
    title = "Porcentaje de mediciones ausentes en cada estación",
    fill = "",
    x = "",
    y = "Estación",
    caption = "MSMK: taller de calidad del aire"
  ) +
  theme(legend.position = "none")

Teniendo en cuenta que una de las variables que consideramos importante en nuestro estudio es pm25, esta variable no tiene una buena calidad en el periodo de los datos anterior al 2005. Por lo tanto, eliminamos aquellas observaciones anteriores al 1 de enero del 2005.

calidad_aire <- calidad_aire %>% 
  filter(fecha >= as.Date("2005-01-01"))

3 Fecha

El conjunto de datos tiene una componente temporal importante. Debemos asegurarnos si están recogidos todos los días del periodo temporal.

Para ello, comenzamos ordenando el conjunto de datos por la variable fecha

calidad_aire <- calidad_aire %>% 
  arrange(fecha)

Y comprobamos que la diferencia entre dos días consecutivos es siempre 1

table(calidad_aire$fecha - lag(calidad_aire$fecha), useNA = "always")
## 
##    1 <NA> 
## 5508    1

Efectivamente, no hay días faltantes en los datos.

4 Resumen de las variables

En la siguienta tabla se muestran estadísticos relevantes para cada variable.

skimr::skim(calidad_aire)
Data summary
Name calidad_aire
Number of rows 5509
Number of columns 18
_______________________
Column type frequency:
Date 1
numeric 17
________________________
Group variables None

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
fecha 0 1 2005-01-01 2020-01-31 2012-07-17 5509

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
anyo 0 1 2012.05 4.35 2005.00 2008.00 2012.00 2016.00 2020.00 ▇▆▆▆▅
mes 0 1 6.49 3.46 1.00 3.00 7.00 10.00 12.00 ▇▅▅▅▇
dia 0 1 15.73 8.80 1.00 8.00 16.00 23.00 31.00 ▇▇▇▇▆
so2 0 1 8.17 4.31 1.60 5.40 7.40 9.56 36.38 ▇▃▁▁▁
co 0 1 0.40 0.18 0.14 0.28 0.34 0.46 1.73 ▇▂▁▁▁
no 0 1 26.72 31.38 1.58 8.08 14.42 31.04 259.69 ▇▁▁▁▁
no2 0 1 42.75 18.34 7.92 29.08 39.77 53.73 133.92 ▆▇▃▁▁
nox 0 1 83.71 64.32 11.50 41.88 62.21 101.42 532.00 ▇▂▁▁▁
pm25 10 1 11.98 6.17 1.50 7.50 10.67 15.00 58.67 ▇▃▁▁▁
pm10 0 1 23.34 13.88 2.58 13.83 20.50 29.17 216.67 ▇▁▁▁▁
o3 0 1 46.23 22.46 2.57 28.00 47.79 63.64 110.07 ▅▆▇▅▁
tol 20 1 3.15 2.43 0.20 1.57 2.50 3.93 31.20 ▇▁▁▁▁
ben 19 1 0.64 0.44 0.12 0.35 0.52 0.78 8.60 ▇▁▁▁▁
ebe 19 1 0.76 0.51 0.10 0.38 0.72 0.95 5.70 ▇▁▁▁▁
tch 0 1 1.42 0.12 1.06 1.33 1.39 1.47 2.42 ▃▇▁▁▁
ch4 0 1 1.25 0.11 0.83 1.18 1.24 1.30 2.16 ▁▇▁▁▁
nmhc 1 1 0.17 0.09 0.01 0.11 0.16 0.21 0.66 ▆▇▁▁▁

En un proyecto real, habría que inspeccionar el comportamiento variable a variable (en este caso es asumible analizar 18 variables). Aquí solamente analizaremos la variable pm25.

4.1 Variable pm25

La distribución de la variable pm25 es la siguiente:

calidad_aire %>% 
  ggplot(
    aes(x = pm25)
  ) +
  geom_histogram() +
  labs(
    title = "Distribución de la variable pm2.5",
    x = "",
    y = "",
    caption = "MSMK"
  ) +
  theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 10 rows containing non-finite values (stat_bin).

La medición de la magnitud a lo largo del tiempo se puede ver en el siguiente gráfico

p <- calidad_aire %>% 
  ggplot(
    aes(x = fecha, 
        y = pm25)
  ) +
  geom_line() +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y")
p

La información del gráfico anterior es difícil de asimilar. En este caso podría ayudarnos si tuviésemos un gráfico interactivo. Podemos convertir un gráfico de ggplot2 a interactivo con ayuda del paquete plotly

plotly::ggplotly(p)

En el gráfico anterior se aprecia vagamente un comportamiento cíclico de forma que los valores parecen incrementarse en los meses centrales del año. Podemos analizar este comportamiento en el siguiente gráfico

calidad_aire %>% 
  ggplot(aes(x = mes, y = pm25)) + 
  geom_jitter(alpha = 0.25) +
  geom_boxplot(aes(group = mes), alpha = 0, color = "firebrick", size = .75) +
  scale_x_continuous(breaks = 1:12)
## Warning: Removed 10 rows containing non-finite values (stat_boxplot).
## Warning: Removed 10 rows containing missing values (geom_point).

A priori, el mes parece relevante en el comportamiento de la variable. En la fase de modelización comprobaremos si influye en la predicción o no.

5 Preparación para modelización

Antes de modelizar es necesario abordar dos tareas:

  1. Generar nuevas variables a partir de la información ya existente en los datos o de información externa.
  2. Modificar las variables para que el modelo solo pueda utilizar la información que tendría en el entorno real de predicción.

5.1 Generación de variables

El conocimiento que se tenga del contexto del proyecto en muchos casos es fundamental para obtener mejores modelos. Distinguimos dos tipos de generación de variables. El primero de ellos es la generación variables a partir de la información que contiene el propio conjunto de datos. El segundo será introducir información externa.

5.1.1 Información del propio conjunto de datos

5.1.1.1 Día de la semana (dia_sem y finde)

Al tratarse de datos de calidad del aire, puede ser un factor relevante el día de la semana. A priori (y será algo que tengamos que comprobar en la fase de modelización), el mayor o menor uso de los vehículos, podría influir. Generamos la nueva variable dia_sem.

calidad_aire$dia_sem <- 
  lubridate::wday(calidad_aire$fecha, 
                  label = TRUE,
                  week_start = 1)

Igual que hicimos anteriormente con el mes, podemos representar el comportamiento de la variables pm25 en función del día de la semana.

calidad_aire %>% 
  ggplot(aes(x = dia_sem, y = pm25)) + 
  geom_jitter(alpha = 0.25) +
  geom_boxplot(aes(group = dia_sem), alpha = 0, color = "firebrick", size = .75)
## Warning: Removed 10 rows containing non-finite values (stat_boxplot).
## Warning: Removed 10 rows containing missing values (geom_point).

En el gráfico anterior, el sábado y el domingo están, en media, ligeramente por debajo de los otros días. Por este motivo, vamos a crear también una variable que indique si la observación se corresponde con el fin de semana.

calidad_aire$finde <- as.numeric(calidad_aire$dia_sem %in% c("sáb", "dom"))

5.1.1.2 Valor medio mensual (*_med_mes)

Aunque disponemos de la variable mes y el modelo de predicción la puede tener en cuenta, a veces ayuda a tener el valor medio de cada mes

calidad_aire <- calidad_aire %>% 
  group_by(mes) %>% 
  mutate_at(vars(so2:nmhc), list(med_mes = mean), na.rm = TRUE) %>% 
  ungroup()

5.1.1.3 Incremento de mediciones (*_inc1)

Es habitual en el caso de series temporales incluir el incremento porcentual respecto del día anterior.

calidad_aire <- calidad_aire %>% 
  mutate_at(vars(so2:nmhc),
            list(inc1 = function(x) (x - lag(x))/lag(x))
            )

Se podrían generar muchas más variables y, conforme se va ganando experiencia en un determinado contexto, mejor sabremos qué variables suelen ser relevantes.

5.2 Estructura adecuada para la modelización

Recordemos que el objetivo es predecir el valor de pm25 en el día posterior. Eso significa que, si queremos predecir el valor para el 9 de febrero, tenemos que suponer que solamente conocemos la información hasta el día 8 de febrero. Por lo tanto, el valor que podremos utilizar de las variables para predecir un día, debe ser la información disponible hasta el día anterior. Por ejemplo, no podemos utilizar el valor de so2 en el mismo día que el valor de pm25 que queremos predecir. Necesitamos que el valor de cada variable esté retrasado en un día. Esto lo podemos hacer mediante la función lag. Para entenderlo, vamos hacer primero un ejemplo. Si tuviésemos el vector c(1,2,3,4), si aplicamos la función lag obtendríamos

lag(c(1,2,3,4))
## [1] NA  1  2  3

Obviamente, el valor anterior del primer elemento del vector es desconocido y por eso aparece como NA.

Ahora generamos una nueva variable _lag por cada variable que tengamos que retrasar. Por ejemplo

calidad_aire$so2_lag <- lag(calidad_aire$so2)

Este procedimiento habría que repetirlo demasiadas veces y sería demasiado pesado para hacerlo de forma manual. Para automatizarlo podemos utilizar la función mutate_at:

calidad_aire <- calidad_aire %>% 
  mutate_at(vars(so2:nmhc, so2_med_mes:nmhc_inc1), list(lag = lag))

Nota: lo que acabamos de hacer se puede traducir como: aplica la función lag a aquellas variables que están entre so2 y nmhc. Al utilizar list(lag = lag), cada variable que se crea termina en _lag.

Para finalizar la creación de estas variables, debemos eliminar todas las variables que no son lag y quedarnos solamente con pm25 que es la variable que queremos predecir.

calidad_aire <- calidad_aire %>% 
  select(fecha, anyo:dia, dia_sem, finde, ends_with("_lag"), pm25)

Es interesante conocer la correlación de la variable objetivo con respecto a las predictoras:

cor_pm25 <- cor(calidad_aire$pm25, 
    select(calidad_aire, ends_with("lag")), 
    use = "complete.obs"
    )

sort(cor_pm25[1,])
##           o3_lag      o3_inc1_lag   o3_med_mes_lag  so2_med_mes_lag 
##      -0.26097887      -0.15548085      -0.09729124       0.09123486 
##   co_med_mes_lag pm10_med_mes_lag  no2_med_mes_lag nmhc_med_mes_lag 
##       0.10186458       0.10445794       0.11884661       0.11903891 
##     ebe_inc1_lag  nox_med_mes_lag   no_med_mes_lag  ben_med_mes_lag 
##       0.12204174       0.12592394       0.12718846       0.12758728 
##  ch4_med_mes_lag  tch_med_mes_lag  tol_med_mes_lag  ebe_med_mes_lag 
##       0.12866286       0.12872078       0.12938145       0.13110331 
##     ch4_inc1_lag    nmhc_inc1_lag      no_inc1_lag     tol_inc1_lag 
##       0.14396524       0.15538853       0.16101512       0.16662015 
##     no2_inc1_lag          ch4_lag     tch_inc1_lag     so2_inc1_lag 
##       0.17772042       0.17826189       0.18764322       0.19016680 
##     nox_inc1_lag pm25_med_mes_lag         nmhc_lag     ben_inc1_lag 
##       0.19828819       0.20517579       0.20549780       0.22828950 
##      co_inc1_lag    pm25_inc1_lag          tch_lag    pm10_inc1_lag 
##       0.25573605       0.28546376       0.29525896       0.30514457 
##          ebe_lag          so2_lag          ben_lag          tol_lag 
##       0.36193726       0.39123372       0.39895123       0.48218583 
##           no_lag           co_lag          nox_lag          no2_lag 
##       0.49360402       0.51821069       0.52853479       0.55917953 
##         pm10_lag         pm25_lag 
##       0.61601447       0.69459147

Nota: en la correlación utilizamos use = "complete.obs" para que no tenga en cuenta los NA en el cálculo.

Puedes ver que la mayor correlación de pm25 se da con pm25_lag. Podemos corroborarlo representando la relación entre ambas variables mediante un diagrama de dispersión:

calidad_aire %>% 
  ggplot(
    aes(x = pm25_lag,
        y = pm25
    )
  ) +
  geom_point(alpha = 0.25) +
  geom_smooth() +
  theme_minimal() +
  labs(
    title = "Relación de PM2.5 con el día anterior",
    x = "PM25 del día anterior",
    y = "",
    caption = "MSMK: taller de calidad del aire"
  )
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 15 rows containing non-finite values (stat_smooth).
## Warning: Removed 15 rows containing missing values (geom_point).

6 Train y test

Por último, como hacemos habitualmente, vamos a dividir el conjunto de datos en train y test. Entrenaremos con datos hasta 2019-09-01 y los restantes para test:

train <- calidad_aire[calidad_aire$fecha < as.Date("2019-09-01"),]
test <- calidad_aire[calidad_aire$fecha >= as.Date("2019-09-01"),]

7 Exportación de la información

Igual que hicimos en la fase anterior, almacenamos estos datos en el disco duro.

saveRDS(train, file = "data/train.RDS")
saveRDS(test, file = "data/test.RDS")